arXiv: 1501.02005vl [cond-mat.str-el] 8 Jan 2015 


Propagation of a single hole defect in the one-dimensional Bose-Hubbard model 
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We study nonequilibrium dynamics in the one-dimensional Bose-Hubbard model starting from an 
initial product state with one boson per site and a single hole defect. We find that for parameters 
close to the quantum critical point, the hole splits into a core showing a very slow diffusive dynamics, 
and a fast mode which propagates ballistically. Using an effective fermionic model at large Hubbard 
interactions U, we show that the ballistic mode is a consequence of an interference between slow 
holon and fast doublon dynamics, which occurs once the hole defect starts propagating into the 
bosonic background at unit filling. Within this model, the signal velocity of the fast ballistic mode 
is given by the maximum slope of the dispersion of the doublon quasiparticle in good agreement 
with the numerical data. Furthermore, we contrast this global quench with the dynamics of a single 
hole defect in the ground state of the Bose-Hubbard model and show that the dynamics in the latter 
case is very different even for large values of U. This can also be seen in the entanglement entropy 
of the time evolved states which grows much more rapidly in time in the global quench case. 

PACS numbers: 05.30.Jp, 05.70.Ln, 67.85.-d 


I. INTRODUCTION 

An investigation of the states of matter and excita¬ 
tions of a quantum system is usually based on its dynam¬ 
ical response to external perturbations. This includes 
frequency-dependent response functions and scattering 
cross sections tested in condensed matter experiments, 1,2 
as well as time-dependent correlation functions measured 
in experiments on cold atomic gases. 3-10 The investiga¬ 
tion of time-dependent correlation functions in and out of 
equilibrium has recently attracted renewed interest due 
to experiments on cold atomic gases in optical lattices 
which provide useful schemes for the preparation of quan¬ 
tum states and the tuning of interactions. 10 ” 12 Recent 
developments allow, in particular, for single-site address¬ 
ability, both in preparation and detection. 13 ” 1 ' 

A lot of attention has been focussed on the dynamics 
of quasi one-dimensional systems where quantum fluctu¬ 
ations are strong and can prevent the ground state from 
attaining any type of long-range order. In this case, the 
system can often be described by Luttinger liquid (LL) 
theory, 18 a hydrodynamic approach for the collective ex¬ 
citations (density fluctuations), which propagate with an 
interaction-dependent velocity. Recent studies, however, 
have shown that LL theory is not sufficient to fully de¬ 
scribe the dynamics of such systems even at zero temper¬ 
ature. Nonlinearities of the dispersion have to be taken 
into account, which lead to x-ray edge-type singularities 
in response functions and new quasiparticles. An exten¬ 
sion of standard LL theory, the so-called nonlinear Lut¬ 
tinger liquid (NLL) theory, has been developed, which is 
capable of describing the collective excitations and the 
new quasiparticles on an equal footing. 19 This allows, 
for example, to calculate the dynamical structure fac¬ 
tor for the Lieb-Liniger model 20 —the continuum limit 
of the Bose-Hubbard model in the superfluid phase— 


which shows power-law singularities at the edges of its 
support. 19 Very recently, it has been demonstrated that 
the NLL approach can also be extended to small finite 
temperatures. 21 One of the main findings of the latter 
study is that the combination of quantum and ther¬ 
mal fluctuations can dramatically alter the dynamical 
response compared to the zero temperature case. For 
the spin-1/2 XXZ chain, in particular, the autocorrela¬ 
tion function at finite temperature and long times was 
shown to exhibit a slow diffusive decay ~ I/'/i in time 
t, explaining the diffusive response measured in NMR 
and /rSR experiments. 22 ” 24 To also understand spin and 
energy transport in the linear response regime, one has 
to identify, in addition, all the local and quasilocal con¬ 
served charges of the microscopic model under considera¬ 
tion. For the XXZ chain, for example, the spin current is 
only partly conserved, so that a ballistic contribution co¬ 
exists with a diffusive contribution. 22,23 This result shows 
that even the linear response in a very simple strongly 
correlated one-dimensional quantum model can already 
be quite complicated and very different from phenomeno¬ 
logical theories. 25 

Dynamics and transport far from equilibrium are, on 
the other hand, much less understood so far even for 
the simplest interacting models. Much of the interest 
has concentrated on quantum quenches, i. e., the dy¬ 
namics ensuing after preparing the system in the ground 
state of a Hamiltonian, and then suddenly changing one 
of its microscopic parameters. The expectation values 
of local observables at long times after the quench in 
a generic, infinitely large quantum system are expected 
to become time independent and to be described by a 
Gibbs ensemble. 26 For integrable quantum systems, on 
the other hand, the additional local and quasilocal con¬ 
served charges have to be taken into account as well, 
leading to a generalized Gibbs ensemble. 27 ” 29 While (gen- 



2 


eralized) Gibbs ensembles can help us to calculate the 
equilibrated values of local observables after a quantum 
quench for a system in the thermodynamic limit, they 
are not helpful in identifying the dominant relaxation 
processes and in describing the nonequilibrium dynamics 
itself. A notable known feature in the quench dynamics 
of a model with short-range interactions, starting from 
an initial state with a finite correlation length, is the 
existence of Lieb-Robinson bounds which state that cor¬ 
relations can only spread with a finite speed. 30,31 The 
emergence of an effective light cone has also been ob¬ 
served in experiments. 5,32 Furthermore, such bounds can 
exist for finite-temperature initial states as well. 33 

A possible way to systematically probe nonequilibrium 
dynamics and transport beyond linear response is to con¬ 
sider local perturbations. 6,7,34-38 The advantage of this 
approach is that local perturbations can easily be cre¬ 
ated in an experimental setup of cold atoms, and that 
one can follow their evolution in time by measuring one- 
point functions instead of two-point correlators. For the 
Bose-Hubbard model, for example, the time evolution of 
local Gaussian density perturbations on top of the ground 
state has been studied by numerical density matrix renor¬ 
malization group (DMRG) methods. 34 In this work, it 
has been shown that for weak and intermediate inter¬ 
actions, small density perturbations travel through the 
system at a velocity which is consistent with the sound 
velocity in the corresponding continuum, exactly solv¬ 
able Lieb-Liniger model. The expansion of density-wave 
packets in the one-dimensional Bose-Hubbard model has 
also been studied experimentally in an out-of-equilibrium 
situation. 39 Here, a gas of 39 K atoms was prepared in a 
harmonic confinement. The confinement was then low¬ 
ered so that the atoms started to expand into the empty 
optical lattice. A slowing down of the expansion of the 
cloud was observed near the critical point in the ground 
state phase diagram. The same experiment in two dimen¬ 
sions showed, in addition, a separation of the cloud into 
a diffusive core and a background which spreads ballisti- 
cally. An important feature of such expansion dynamics 
is that more than one velocity can be present due to the 
formation of bound states. Examples for such behavior 
are the XXZ model at large repulsive nearest-neighbor 
interaction, 6,40-42 as well as the Bose-Hubbard model, 
where the onsite repulsive interaction can lead to effec¬ 
tively bound states of doublons. 43,44 

In this paper, we will study the non-equilibrium dy¬ 
namics of one of the simplest local perturbations: a sin¬ 
gle hole defect in an initial product state with one bo¬ 
son per lattice site under the unitary time evolution of a 
Bose-Hubbard Hamiltonian. This setup can easily be re¬ 
alized using cold atomic gases, which offer single-site ad¬ 
dressability in the preparation of initial states and in the 
detection of bosons after a quantum quench. 13-17 Here, 
the hole defect spreads into a background which itself 
displays nontrivial dynamics, contrasting our setup from 
the case of an expansion into an empty lattice, which 
has been studied thoroughly for the bosonic 39,43-48 and 


fermionic 49-53 versions of the Hubbard model. We will 
show, in particular, that the defect will either spread bal- 
listically through the lattice or separate into a diffusive 
core coexisting with a fast ballistic mode, depending on 
the strength of the Hubbard interaction U. This phe¬ 
nomenon is reminiscent of the coexistence of ballistic and 
diffusive transport in the XXZ model at finite tempera¬ 
ture, however, here it is not related to local or quasilocal 
conserved charges, but is rather, as we will show, a con¬ 
sequence of an interference effect between slow holon and 
fast doublon quasiparticle dynamics. 

Our paper is organized as follows: In Sec. II we de¬ 
fine the model, the initial states, and the time-dependent 
correlation functions which we study in the following. In 
Sec. Ill we explain the numerical method used to simulate 
the quench dynamics and introduce an effective fermionic 
model 54 to approximate the dynamics for large interac¬ 
tions. In Sec. IV we present the results of the numerical 
calculations and analyze these results using the effective 
fermionic model. In Sec. V we summarize our main re¬ 
sults and discuss some of the remaining open questions. 


II. MODEL 

We consider the one-dimensional Bose-Hubbard 
model, 

H = J Y^( b l b j +1 + hx -) + \ ^l n A n i - !)> (!) 

3 3 

where bj (bj) creates (annihilates) a boson at site j, 

rij = b^bj is the local density operator, U is the strength 
of the on-site repulsive interaction, and in the following, 
we set J = 1. At unit filling, this model exhibits a quan¬ 
tum phase transition at U c ss 3.37 between a superfluid 
and a Mott insulating state. 55 We are interested in the 
time-dependent density profile C r (j,t) after removing a 
single particle at the center of the chain from a prepared 
uniform initial state |\E' r ) 

C r (j, t) = ('b r \bj— 0 e iHt rije~ iHt bj—o\'i’ r ). (2) 

Our main focus is on the time evolution starting from the 
initial product state v I / r = JJ ) = J]^6)|0}. To understand 
the unique features of the non-equilibrium dynamics, we 
will compare our results with the equilibrium correlation 
function obtained by using the ground state |'F, r = s ) of 
the time-evolving Hamiltonian H at unit filling as the 
initial state in Eq. (2). Alternatively, the latter case can 
be viewed as a local quench—defined as a local modifica¬ 
tion of the Hamiltonian or the initial state—while in the 
former case, an additional global quench from the ground 
state at infinite repulsion is performed. 
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III. METHOD 
A. Numerics 

We calculate the correlation function (2) using the light 
cone renormalization group (LCRG) algorithm. 56 This 
algorithm is based on the time-dependent DMRG 57 ~ 60 
and yields results directly in the thermodynamic limit. 
In our simulations we restrict the maximum number of 
bosons per lattice site to 5 for U < 6, and to 3 for 
U > 6. By comparing with simulations where the max¬ 
imum number of bosons is smaller we make sure that 
the results are converged on the scale of the figures we 
present in the following. Furthermore, we make sure that 
the discarded weight in each renormalization group step 
is smaller than ICE 8 . For the simulation times shown, 
this requires to increase the number of states kept in the 
truncated Hilbert space up to 8 000. 



j 

FIG. 1: (Color online) The density profile C p (j,t) = 1 — 
J^(2t) for U = 0 and U = oo. The dashed lines show the 
linear spreading of the defect density, j = ivt, with a defect 
velocity v = 2. 

For the calculations where the initial state is the 
ground state |'F g ), we first perform an imaginary time 
evolution with the correct chemical potential ^ to project 
an arbitrary starting vector onto the ground state at 
unit filling. The chemical potential is required because a 
canonical simulation in the thermodynamic limit is im¬ 
possible. We check that |(’F g |n J j’I' £ ,) — 1| < 10 -6 , i. e. we 
ensure convergence in the density to more than 6 digits. 

B. Analytical results and approximations 

The Bose-Hubbard model, Eq. (1), has two integrable 
limits: the limit of free bosons, [7 = 0, and the limit of 
hard-core bosons, U = oo. In both cases, we are dealing 
with single-particle dynamics and the density profile can 
be calculated exactly: 

CZ=°(j,t) = Cp =oc (j, t) = 1 - 


where Jj is the jth Bessel function of the first kind. The 
density profile for these cases is shown in Fig. 1: The de¬ 
fect density spreads in a light cone with a velocity v = 2. 
If we create a doublon instead of a hole defect, i.e., in¬ 
terchange bj =0 bj =o hi Eq. (2), the corresponding cor¬ 
relation function C p is given by C p =0 (j,t) = 1 + 

and C p =oc (j,t) = 1 + In this case, the velocity 

of the expansion is given by v = 2 for U = 0 and v = 4 
for U = oo. The by a factor of 2 larger velocity in the 
latter case is a consequence of the Bose factor. 

In the following, we will call any dynamics similar to 
the one shown in Fig. 1 ballistic. For a more precise 
definition, we follow Refs. 36,48 and consider 

OO 

R\t) = E (4) 

j=—o o 

where R(t) is a weighted measure for the distance over 
which the defect density has spread at time t. For the 
free and hardcore boson cases, the sum can be calculated 
exactly, using Eq. (3) and a standard identity for Bessel 
functions, resulting in 

OO 

R 2 (t)= E fJh‘*t) = 2t 2 . (5) 

j = - oo 

The functional dependence of R 2 on time is thus the same 
as for a linearly dispersing, ballistic particle, 

oo *2 

R 2 (t) = y [S(j + vt) + S(j - vt)] =v 2 t 2 . (6) 

j =—oo 

For large but finite U, an exact solution is not possible. 
Since high occupancies are suppressed, we can, however, 
approximate the dynamics in this limit by restricting the 
local Hilbert space to three states only: |0)j, l)j, and 
12 )j. We can now interpret the state |l)j as the vacuum, 
the doublon state |2 )j as a spin-up particle, and the holon 
state |0)j as a spin-down particle. This is made explicit 
by the mapping 

bj = \[2btf + (7) 

where bj a are hardcore bosons. 54 The hardcore constraint 
for identical particles can be fulfilled by a second mapping 
onto fermions 

= ^j c jt > bji = Zje ( 8 ) 

where Zj = J}., . nja is a Jordan-Wigner string, 

Cj a are fermionic operators, and rij a = c { , jCT Cj a . Applying 
both transformations maps the bosonic Hamiltonian (1) 
in the restricted Hilbert space onto the fermionic Hamil¬ 
tonian in momentum space, 

H = Y^ Ek ° C kv C *° + i ^ A k( c U c ll.- C nCkt) 

key k 



k,k' ,q 


( 3 ) 


( 9 ) 
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with V —> oo. Here, we have defined k = —k, A k = 
2\/2sinfc, and E ka = U/2 — (3 + cr) cos k, where cr = ±1 
for the spin-up and spin-down dispersions, respectively. 
The infinite potential V projects out the unphysical state 
|j'j,)j, i.e., the state in which a doublon and a holon are on 
the same site. Within the truncated local Hilbert space, 
this transformation is exact. 

As a further approximation, we will ignore the con¬ 
straint and allow for double occupancies by setting V = 
0, leading to the unconstrained fermion (UF) model con¬ 
sidered in Ref. 54. This is justified as long as the density 
of doublons and holons is low, which is, at large enough 
interaction U, the case both for the ground state and the 
dynamics starting from the product state. The Hamilto¬ 
nian, Eq. (9), can then be diagonalized by a Bogoliubov 
transformation, 

Cka = Ika COS 0 k + i%- sin 9 k , (10) 

with Bogoliubov angle 9 k = arctan[2Afc/(£’^ + E’^)]/2. 
This leads to a diagonal Hamiltonian in the Bogoliubov 
quasiparticles, H = ]U fccr e k <rl kcr lkai with dispersions 

£fccr = - {Ek\ + E k |) 2 + 4A| + — (E k f — E k ±). (11) 

The ground state |of the Bose-Hubbard model 
( 1 ) is then approximated by the vacuum of the 7 - 
quasiparticles, 7 fc (J |'I' ff ) = 0 , while the product state |^ p ) 
takes the form of a BCS state, 

\^ p ) = II (cos^fe +* sin6 'fc7l t 7| 4 .) I^ 9 )- (12) 

k 

Note that the doublon quasiparticle dispersion becomes 
gapless at XJ = 8 , eof = 0, showing that the UF approx¬ 
imation breaks down for smaller interactions. Instead 
of setting V = 0 in Eq. (9), the constraint can be taken 


into account approximately by a summation of ladder di¬ 
agrams which leads to a Bethe-Salpeter equation . 61 The 
corresponding self-energies will then shift both quasipar¬ 
ticle dispersions e ka up, so that the model remains stable 
even below U = 8 . However, this formalism cannot di¬ 
rectly be used to calculate the nonequilibrium correlation 
function we are interested in. We will therefore only use 
the UF model to explain qualitative features of the defect 
dynamics for interactions U > 8 . 


IV. RESULTS 
A. Large interactions 

In Fig. 2 (a) and (b), we show the numerically obtained 
profiles C p (j,t) for large, but finite values of the interac¬ 
tion strength. In both cases, a clear qualitative change as 
compared to the U = 00 case shown in Fig. 1 is observed. 
Instead of a ballistic spreading with a single velocity, we 
find a separation into a central light cone spreading with 
velocity v p , and an extra contribution which moves with 
a faster velocity v p and clearly separates itself from the 
main defect density. To extract the Lieb-Robinson ve¬ 
locity of the first signal, i>®, from the numerical data, we 
take the point j(t) where the defect density has reached 
half of the height of the first peak of its time-dependent 
distribution as reference point, and perform a linear fit 
according to j(t) = v^t. 56 

Next, we show that the presence of two velocities in 
the defect dynamics can be explained using the effec¬ 
tive UF model. Expressing the bosonic density operator 
in the restricted Hilbert space in terms of the fermions, 
rij = 1 + Cj^Cj-f — Cj^Cji, the correlation function C p can 
be calculated within the UF model approximation, see 
Appendix A for details. We obtain 


C P {j,t) = 1 - 


■j^ /* 7T 2 -j^ /* 7T 

— J dk e lk3 (g k i(t) cos 2 9 k + g k ^(t) sin 2 9 k ) - — J dke rk ° cos9 k sm.9 k (g kir (t) - g k f(t)) 


1 - Jf(2t) ( 1 - 


2 j (j - 2 ) 
UH 2 


U 2 


- jjm 4t) + - ~ t j j+ i(4t) 


j(j + 1 ) 


(13) 


with propagators g ka (t) = bfc^bL) 0 )) = e lekat . and 
we have performed the continuum limit. In the second 
line the integrals have been evaluated within an expan¬ 
sion in l/U up to order U~ 2 . The result shows that 
contributions propagating with two different velocities as 
well as interference terms are present in the dynamics. 
The two different velocities can be determined from the 
maximal slope of the two dispersion relations, 


dt k (7 

v ° = m r -w (14) 

From the l/t/-expansion in Eq. (13) we see that 17 —> 4 
and 17 —► 2 for U —> 00 . This is consistent with the 
discussion of the holon and doublon dynamics at U = 00 
inSec. Ill B. 

While the density profiles for the UF model, Fig. 2 
(c) and (d), quantitatively differ from the LCRG results, 
they clearly show the same qualitative features. In partic- 
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FIG. 2: (color online) Density profiles C p / g (j,t) for large U: (a,b) Numerical results by LCRG for the initial product state 
|flip). (c,d) Results for the product state by numerically evaluating the integrals (13) in the UF model approximation. (e,f) 
LCRG results with the ground state |'F 9 ) as initial state. The dashed lines in (a,b) and (e,f) denote the numerically extracted 
velocities v p and v a g , respectively. The dotted lines in (a,b) are a guide to the eye indicating the slow core velocity v p . The 
dashed lines in (c,d) denote the velocities v a defined in Eq. (14). 


ular, we find good agreement between the renormalized 
holon velocity and the slow core velocity, v±, ~ v p , and 
the renormalized doublon velocity clearly determines the 
velocity of the first signal, u-j- w v p . In Fig. 3 we show 
the contributions to C p , Eq. (13), which contain the dou¬ 
blon propagator separately. Comparing with Fig. 2 
(a) and (b) we see that these contributions are indeed 
responsible for the fast mode. From the analytical re¬ 
sult in the 1/[/-expansion, see second line in Eq. (13), we 
can infer the weight W of these contributions to leading 
order, 


w ^ ST' 2 j(j — 2 ) 2 4 

lJ2 t 2 J j( 2t >=lJ2- ( 15 ) 

j—~ oo 

Another way to see that the mixing of holon and dou¬ 
blon excitations is responsible for the fast mode is to 
compare with the ground state correlation function C g , 
shown in Fig. 2 (e) and (f), where the fast mode is absent, 
and only a single velocity v g can be observed. Within the 
UF model the initial state is now the vacuum of the 7 - 
quasiparticles. From Eqs. ( 8 ) and (10) one immediately 
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FIG. 3: (color online) Contributions to C p (j, t) containing the 
gk-f propagator, see Eq. (13), for (a) U = 10 and (b) [7 = 8. 


sees that in this case the < 7 /--(- propagator does not show 
up in the correlation function C g \ accordingly, we see 
that v g ~ 14 . The density profile therefore evolves dif¬ 
ferently in time, as is exemplarily shown in Fig. 4, where 
cuts of the density profiles for fixed times for the corre¬ 
lation functions C p and C g are compared. The difference 
in the defect dynamics starting from the ground state as 
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opposed to the product state is also clearly seen in the 
entanglement entropy, S en t(t), between two halves of the 
infinite chain where the defect initially sits at the bound¬ 
ary, see Fig. 5. The extraction of a particle from the 
ground state of the system constitutes a local quench. In 
this case the entanglement entropy grows very slowly in 
time and the numerical simulation can be performed for 
long times . 62 For the initial product state, however, in ad¬ 
dition to the local quench, a global quench is performed, 
because the initial product state I’Fp) only becomes the 
exact ground state at U = oo. Hence, we see the well- 
known behavior for a global quench of the entanglement 
entropy increasing quickly and linearly with time . 30,63 



FIG. 4: (color online) LCRG results for U = 8 : (a) C p (j,t) 
compared to (b) C g (j,t). In (a) the position of the fast mode 
is marked by arrows. Subsequent curves are shifted by 0.1. 



FIG. 5: (color online) Entanglement entropy S ent (t) for U = 
10 between two halves of the infinite chain where the defect 
initially sits at the boundary. The solid line shows the linear 
increase of S en t(t) for the defect in the initial product state 
|$ p ), as expected for a global quench, while the dashed line 
represents the much slower increase of S en t(t) for the defect 
in the initial ground state |4/ 9 ) (local quench). 


B. Intermediate and weak interaction 

For interaction strengths U < 8 we can no longer use 
the UF model and have to fully rely on numerical results. 
From the numerical data shown in Fig. 6 (a-d), we see 
that the separation between the fast mode and the core 


becomes even more pronounced for intermediate interac¬ 
tion strengths. Particularly interesting is the case U = 3, 
shown in Fig. 6 (c,d), where the core is almost frozen in 
and separated from the fast mode by a region where the 
density is even larger than one. While the fast mode still 
spreads ballistically, we can no longer assign a velocity 
to the core region, hinting at diffusive transport. The 
dynamics is thus very different from that of a defect in 
the ground state, shown in Fig. 7, which remains purely 
ballistic. For even smaller U, shown in Fig. 6 (e,f), the 
signal velocity u® decreases further so that a separation 
from a core region can no longer be resolved, at least not 
on the numerically accessible time scales. Eventually, 
purely ballistic dynamics is restored at U = 0, where the 
density profile is again given by Eq. (3), see Fig. 1. 

To quantify this change in the expansion dynamics of 
the defect for the product state |T p ), we calculate R 2 {t), 
as defined in Eq. (4). The results for various values of U 
are shown in Fig. 8 (a). At short times t < 1 we see that 
the dynamics is almost independent of the interaction 
strength and the defect spreads ballistically, R 2 ~ t 2 . 
At longer times there is a clear difference between large 
and intermediate U values. While for (7 = 9 the dy¬ 
namics is also ballistic at intermediate times according 
to this criterion, the curves for U = 4.5 and U = 3 be¬ 
come linear, R 2 ~ Dt, indicating a dominantly diffusive 
spreading of the defect at this time scale with a diffusion 
constant D. 36,48 At these intermediate timescales we can 
extract the diffusion constant D by linear fits; the results 
are shown in Fig. 8 (b). We see that the spreading of 
the core depends on interaction strength U and is slow¬ 
est for values close to the critical point U c « 3.37 which 
separates the superfluid from the Mott insulating regime 
in the ground state phase diagram. A similar slowing 
down near a critical point was also observed in the free 
expansion dynamics of a product state into an empty 
lattice , 39,48 as well as in the non-equilibrium dynamics 
of the XXZ chain . 64,65 Note, however, that the coexis¬ 
tence of diffusive and ballistic transport, which is clearly 
seen in Fig. 6 (a)-(d), cannot be resolved by the simple 
quantity R 2 (t) on intermediate timescales. In addition, 
R 2 (t) will always be dominated by a ballistic contribution 
in the long-time limit, no matter how small this contri¬ 
bution is. 

A comparison of the different velocities observed in the 
dynamics of a single defect is shown in Fig. 8 (c). In con¬ 
trast to the diffusive core, critical slowing down is not 
observed in the velocity of the fast ballistic mode, u®. In¬ 
stead, the velocity v® is a monotonously increasing func¬ 
tion in U with u® = 2 at U = 0, eventually approaching 
the value u® —> 4 for large U as expected for a single 
doublon in the initial product state, see Sec. IIIB. For 
U > 8 , Vp is well approximated by the doublon quasipar¬ 
ticle velocity v-\ determined from the UF model. On the 
other hand, the velocity u® for the first front of the cor¬ 
relation function C g is not monotonous: Here we have 
Vg(U = oo) = 2, and the velocity first increases with 
decreasing U and approximately agrees with the holon 
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FIG. 6: (color online) Density profiles C p (j,t) at weak and intermediate interactions for (a,b) U = 4.5, (c,d) U = 3.0, and (e,f) 
U = 1.5. The dashed lines in (a,c,e) are the signal velocities v p as extracted from the numerical data. In (b,d,f) subsequent 
curves are shifted by 0.1. 
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FIG. 7: (Color online) Density profiles C g (j,t) starting from 
the ground state l^g) for (a) U = 4.5, and (b) [7 = 3. The 
dashed lines indicate the extracted signal velocities v s g . 


Ull 



(16) 


One has to keep in mind, however, that removing a par¬ 
ticle from the ground state is not a low energy excitation 
of the system, but rather involves all momentum modes, 
bj =o = IV -1 ^2 k bk- Only for sufficiently large U in the 
superfluid regime, where the linearization of the spec¬ 
trum is a good approximation for a large region around 
k = 0, can we expect that both velocities are similar. 
This is indeed the case, as seen in Fig. 8(c). We find, fur¬ 
thermore, that the two velocities are also similar for U 
values even slightly above the Mott transition. Here the 
excitation gap is relatively small so that the linear ap¬ 
proximation of the dispersion still works approximately 
for a large range of momenta k > 0. 


quasiparticle velocity v\,, see Eq. (14), in the regime 
where the UF approximation is valid. For smaller values 
of U, Vg decreases with decreasing U. In the superfluid 
regime, we can then compare v s g with the sound velocity 
of the corresponding continuum Lieb-Liniger model, 66,67 


V. SUMMARY AND CONCLUSIONS 

We have considered one of the simplest setups to sys¬ 
tematically study the non-equilibrium dynamics of the 
one-dimensional Bose-Hubbard model: a single hole de- 


























































FIG. 8: (color online) (a) R 2 (t ) for the spreading of the defect 
with |\P p ) as initial state (symbols), a quadratic fit for U = 9 
(red line), and linear fits for U = 4.5 and U = 3 (blue lines), 
(b) Diffusion constants, R 2 (t ) ~ Dt, in the regime where the 
dynamics is dominantly diffusive at intermediate times. The 
expansion is slowest close to the critical point U c ~ 3.37. (c) 
Velocity v p for the fastest signal in the product state dynamics 
(black circles), compared to the velocity of the signal Vg for 
the ground state dynamics (red diamonds). For U > 8 both 
velocities are well approximated by the UF model velocities 
v-f and V|_, Eq. (14), respectively (black solid lines). The blue 
dashed line denotes the sound velocity vll in the Lieb-Liniger 
model. 

feet in an otherwise uniform state. If this initial state is a 
product state at unit filling, we have shown that the de¬ 
fect separates into a fast ballistic mode and a slower core. 
This behavior, which exists for any finite interaction, can 
be explained by an effective model valid for large values 
of U: doublon quasiparticles give rise to it. Using the 
quantity R 2 {t), which characterizes the extent of the re¬ 
gion over which the defect has spread, we found that the 
core shows signs of a diffusive time evolution at inter¬ 
mediate interactions, leading to a coexistence of ballistic 
and diffusive defect dynamics. While the diffusive core 


shows critical slowing down, i. e. the dynamics is slowest 
close to the quantum critical point, the velocity of the 
fast mode does not. An effective theory which describes 
the coexistence of a diffusive core and the ballistic back¬ 
ground in the dynamics and explains the slowing down 
of the diffusion near the quantum critical point remains 
as an open problem. We showed, furthermore, that the 
diffusive spreading of the defect only occurs when start¬ 
ing from the product state. If one considers instead the 
ground state of the Bose-Hubbard model as the initial 
state, then the dynamics is always purely ballistic, and 
the velocity in the superfluid regime is, for U > 1, well de¬ 
scribed by the sound velocity of the Lieb-Liniger model. 
The proposed quench protocol can easily be implemented 
in an ultracold atom experiment using state-of-the-art 
preparation and detection techniques with single-site res¬ 
olution and offers a systematic way to study the dynam¬ 
ics of local perturbations in the Bose-Hubbard model be¬ 
yond the linear response regime. Importantly, this test of 
the non-equilibrium dynamics only requires to measure 
one-point correlation functions, because translational in¬ 
variance of the initial state is broken by the defect. 
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Appendix A: Density profile from the UF model 

Our objective is to calculate the correlation function 

C p (j,t) = (* p \ble iHt n je - iHt b 0 \* p ) (Al) 

within the approximation of the UF model. Using equa¬ 
tions (7), ( 8 ), (10), and (12), we can express the initial 
state and the bosonic density operator, rij = l + ct^c^q — 
in terms of the 7 -quasiparticles, 


& ° 1i *p> = ^ E ^ n ( cos °p +* sin vMj i °7>. 
1 


^ = ^E ei(fe_fc ' b ' 


k,k' 


cos(0fc - 0fc')(7fc t 7fc'T + 7jfc4.7jL) + i sin {0 k - 0*/)(7fc t 7jL + 7fef7fc't) 


(A2) 

(A3) 


We can now calculate C p straightforwardly by using Wick’s theorem. For example, the contribution of the first 
summand in Eq. (A3) is given by 













9 


_L e i(fc k ' )j cos( 0 fc - 0*0/794. n( co ^, “ isin 0 p 7 p 4 . 7 p t ) 7 j [ t (t) 7 fc -t(t) ]J ( C0S( V + * sin 0^7^ 7 ^ 4.) 7 ^ 4. ) 

q,k,k',q' ' P^Q P'^Q' / 

= -^^sin 2 0 fc - e ,(fe - fe )j sin 0 * sin 0 */ cos( 0 fc - 0 fcQfffct(Qgfc't(O, (A 4 ) 

k k,k' 


where we introduced the propagators g ka (t) = 

( lkcr(t)T k(T (0 )}• Here, we can clearly see that the ap¬ 
pearance of propagators is a consequence of the BCS 
structure of the initial product state, which therefore 
are not present if the initial state is the ground state 
|0 7 ). Calculating the remaining three contributions in 
the same way, we obtain the result given in the first line 
of Eq. (13). 

To approximate the integrals, we can perform an ex¬ 
pansion in 1/U resulting in 


cos 0 k 

sin Ok 

9ka{t) 


1 * • 2 7 

1 - jp sm 


2V2 . , 

—y— sm k - 

g —itU /2^—(3 -\-cr)it cos k 


^sin(2fc), 


(A5) 


Using this expansion leads to the result given in the sec¬ 
ond line of Eq. (13). 
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